The following exercises are modified from Chapter 9 of Geocomputation with R by Rovin Lovelace. Each question lists the total number of points. The breakdown of points can be found at the end of each instruction in parentheses. A general grading rubric can be found on the course website.
Please update “author” to list your first and last names and any collaborators (e.g. Ruth Oliver, Friend1, Friend2)
Due by midnight Saturday 2022-10-08
library(sf)
library(raster)
library(dplyr)
library(spData)
library(spDataLarge)
library(tmap)
library(RColorBrewer)
library(grid)
library(leaflet)
These exercises rely on a new data object based on the world and worldbank_df datasets from the **spData* package.
africa = world |>
filter(continent == "Africa", !is.na(iso_a2)) |>
left_join(worldbank_df, by = "iso_a2") |>
dplyr::select(name, subregion, gdpPercap, HDI, pop_growth) |>
st_transform("+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25")
We will also use the zion and nlcd datasets from the spDataLarge package.
zion = st_read((system.file("vector/zion.gpkg", package = "spDataLarge")))
## Reading layer `zion' from data source
## `/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/spDataLarge/vector/zion.gpkg'
## using driver `GPKG'
## Simple feature collection with 1 feature and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
## Projected CRS: UTM Zone 12, Northern Hemisphere
data(nlcd, package = "spDataLarge")
force(nlcd)
## class : RasterLayer
## dimensions : 1359, 1073, 1458207 (nrow, ncol, ncell)
## resolution : 31.5303, 31.52466 (x, y)
## extent : 301903.3, 335735.4, 4111244, 4154086 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 8 (min, max)
## attributes :
## ID colors levels
## from: 1 #476ba0 Water
## to : 8 #bad8ea Wetlands
Create a map showing the geographic distribution of the Human
Development Index (HDI) across Africa.
- use both base graphics (hint: use plot()) and
tmap) (4)
- name two advantages of each based on the experience (3) - name three
other mapping packages and an advantage of each (3)
# Create a map of the geographic distribution of HDI using plot
plot(africa["HDI"])
# Not sure what the error is that I'm getting
# Create a map of the geographic distribution of HDI using tmap
tm_shape(africa) +
tm_polygons(col = "HDI") +
tm_layout(main.title = "HDI by Countries in Africa",
main.title.size = 1,
bg.color = "tan")
Advantages of plot()
The default legend in plot is really well placed and apparent. It makes it easier for the viewer to know what they are looking at.
The default colors used in plot have high contrast which make it easy to distinguish which countries fit into which category.
Advantages of tmap
Tmap feels a lot like ggplot which many r users are highly familiar with. This makes it nice and intuitive when adding layers and elements to the map.
While the default colors on tmap are not as contrasted as plot(), they show a nice scale ranging for light to dark of the same colors (red and orange here), and they clearly indicate to the viewer where high or low HDI occur in Africa.
Three other mapping packages and when you might use them/their advantages:
Use maps to outline countries and label cities with
points
Use mapview to create maps that are easily
interactive
Use ggplot2 to create basic maps especially with
points made up of lat,long data
Extend the tmap created for question 1 so the legend
has three bins: “high” (HDI above 0.7), “medium” (HDI between 0.55 and
0.7), and “low” (HDI below 0.55). (5)
- change the legend title (5)
- change the class labels (5)
- change the color palette (5)
# Add to the map
hdi_africa <- tm_shape(africa) +
tm_polygons(col = "HDI",
breaks = c(0, 0.55, 0.7, 0.75),
labels = c("Low", "Medium", "High"),
title = "Human Density Index",
palette = brewer.pal(n = 4, name = "Blues"),
contrast = 1.7) +
tm_layout(main.title = "HDI by Countries in Africa",
main.title.size = 0.8,
bg.color = "snow2",
legend.position = c(0.01, 0.01))
hdi_africa
Represent Africa’s subregions on the map. (5)
- change the color palette (5)
- change the legend title (5)
- combine this map with the map from question 2 into a single plot
(5)
# Create a map displaying the subregions of Africa
sub_africa <- tm_shape(africa) +
tm_polygons(col = "subregion",
palette = brewer.pal(n = 5, name = "Spectral"),
alpha = 0.9,
title = "Subregion Name") +
tm_layout(main.title = "Subregions of Africa",
main.title.size = 0.8,
bg.color = "snow2",
legend.position = c(0.01, 0.01))
sub_africa
# Place the two maps side by side in one plot
tmap_arrange(hdi_africa, sub_africa)
Create a land cover map of Zion National Park (5)
- change the default colors to match your perception of land cover
categories (5)
- move the map legend outside of the map to improve readability
(5)
- add a scale bar and north arrow and change the position of both to
improve the maps aesthetics (5)
- add the park boundaries on top of the land cover map (5)
- add an inset of Zion’s location in the context of the state of Utah
(5)
- hint: an object representing Utah can be subset from the
us_states dataset)
# First plot the zion polygon
zi_map <- tm_shape(zion) +
tm_borders(lwd = 2,
col = "black")
# Create a color palette that makes sense for landcover
land_palette <- c("#33A2FF", "#B2B6B9", "#E9E3C8",
"#1D5F2A", "#A9E18E", "#E154CE",
"#58360F", "#265364")
# Add the land cover data
lc_map <- tm_shape(nlcd) +
tm_raster(alpha = 1.0,
palette = land_palette,
title = "Land Type") +
tm_layout(legend.outside = TRUE, # Moving the legend
main.title = "Land Cover in Zion NP",
main.title.size = 1) +
tm_compass(type = "4star", # Adding the compass
position = c("0.02", "0.05"),
size = 4) +
tm_scale_bar(breaks = c(0, 2, 4, 6, 8, 10), # Adding the scale bar
position = c("right", "top"),
bg.color = "white") +
tm_graticules() +
zi_map
# Create the boundaries for the inset map
zion$geom
## Geometry set for 1 feature
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
## Projected CRS: UTM Zone 12, Northern Hemisphere
## POLYGON ((314945.2 4115910, 314803.7 4115949, 3...
box <- st_bbox(c(xmin = 302903.1, xmax = 334735.5,
ymin = 4112244, ymax = 4153087),
crs = st_crs(zion)) %>%
st_as_sfc()
# Double check your map of interest
lc_map
## stars object downsampled to 888 by 1125 cells. See tm_shape manual (argument raster.downsample)
# Make the broader map of Utah
ut_geom <- us_states %>%
filter(NAME == "Utah")
ut_map <- tm_shape(ut_geom) +
tm_polygons() +
tm_shape(zion) +
tm_borders(col = "black") +
tm_shape(box) +
tm_borders(col = "black")
ut_map
# Merge the maps
lc_map
## stars object downsampled to 888 by 1125 cells. See tm_shape manual (argument raster.downsample)
print(ut_map, vp = viewport(0.8, 0.27, width = 0.5, height = 0.5))
# I don't feel like the aesthetics look amazing, but I spent a long time to get them to this point!
Create an interactive map of the world (10)
- include a legend (5)
- change the color palette (5)
- bonus: use leaflet insted of tmap
(2)
# Use tmap to make an interactive map of the world with a legend of the continents
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(world) +
tm_polygons(col = "continent",
palette = "RdYlBu")
# Try using leaflet to make the map
leaf_map <- leaflet(world) %>% addTiles()
leaf_map
factpal <- colorFactor("RdYlBu", world$continent, n = 8)
continents <- unique(world$continent)
leaflet(world) %>%
addPolygons(stroke = FALSE, smoothFactor = 0.2, fillOpacity = 1,
color = ~factpal(continent))
# addLegend(colors = factpal,
# values = continents,
# opacity = 1,
# position = "topright")
# Okay I really can't figure out how to get a legend in here :/
Create THE WORST map! AKA a map that breaks all of the rules of legibility/clarity, yet somehow still passes for a map. We’ll vote on the best worst map (worst worst map?) in class.
# This is the map (plus a few edits) I originally thought we were supposed to make for question #3, and I thought it looked terrible. Luckily, I figured out #3 nd now here is my terrible-looking map :)
tmap_mode("plot")
## tmap mode set to plotting
# Recreate map #2 with a gray scale
gray_hdi <- tm_shape(africa) +
tm_polygons(col = "HDI",
breaks = c(0, 0.55, 0.7, 0.75),
labels = c("Low", "Medium", "High"),
title = "Human Density Index",
palette = brewer.pal(n = 4, name = "Greys"),
contrast = 1.7) +
tm_layout(main.title = "HDI by Countries in Africa",
main.title.size = 1,
bg.color = "snow2")
# Recreate map #3 with a low alpha to prep for stacking
trans_sub <- sub_africa <- tm_shape(africa) +
tm_polygons(col = "subregion",
palette = brewer.pal(n = 5, name = "Spectral"),
alpha = 0.2) +
tm_layout(main.title = "Subregions of Africa",
main.title.size = 1,
bg.color = "navy") +
tm_borders(col = "subregion")
gray_hdi +
trans_sub +
tm_layout(main.title = "HDI in subregions of Africa")
## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).